Skip to content

Conversation

@rflamary
Copy link
Collaborator

@rflamary rflamary commented Jan 27, 2020

I just discovered this important bug in the `ot.emd``function and it feels like opening pandoras's box.

I added in the fisrt commit a test that checks that the constraints are satisfied and we can see in the travis build that the test fails. This comes from the fact that the c++ emd solver is sparse in the sens that if a weight is 0 it will discard the bin and solve a smaller subproblem on the histograms that have non-zero weights. It works very well and is very efficient in practice bug it does not return proper dual variables because it sets the value for the bins with 0 weight to 0 which can lead to constraint violation in practice.

What I will do in this PR:

  • Correct the bug: ot.emd should work by default (return proper dual variables/subgradients)
  • add a postprocessing (as default but and ot.emd parameter) to the dual potential so that they are in a way 'centered' because the invariance to a constant can be a big problem in practice when using theme as subgradients. The pre-processing will not change the objective value but provide some stability (and will always return one of the possible solution).

The post-processing is actually necessary when trying to estimate correct dual variable from the optimal sub-problem if we want to avoid bias in the resulting dual potentials (bnasically everyone at 0 on either the left or right potential).

@rflamary rflamary changed the title [WIP] Bug on dual potential (violation constraint when some weights are 0) [WIP] Bug on dual potential (constraint violation when some weights are 0) Jan 27, 2020
@ncourty ncourty self-requested a review January 28, 2020 22:22
Copy link
Collaborator

@ncourty ncourty left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice bug correction @rflamary. Maybe you can be more explicit in the doc why we need this fix (mostly because of the way EMD is coded in C++). In the longer term, should we plan to integrate this patch rather inside the C++ code ?

def estimate_dual_null_weights(alpha0, beta0, a, b, M):
r"""Estimate feasible values for 0-weighted dual potentials
The feasible values are computed efficiently bjt rather coarsely.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

but

\beta_j = \beta_j -v^b_j \quad \text{ if } b_j=0 \text{ and } v^b_j>0
In the end the dual potential are centred using function
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

potentials
are centered

:ref:`center_ot_dual`.
Note that all those updates do not change the objective value of the
solution but provide dual potential that do not violate the constraints.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

potentials

beta0 : (nt,) numpy.ndarray, float64
Target dual potential
a : (ns,) numpy.ndarray, float64
Source histogram (uniform weight if empty list)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

histogram -> distribution ?
weight -> weights

a : (ns,) numpy.ndarray, float64
Source histogram (uniform weight if empty list)
b : (nt,) numpy.ndarray, float64
Target histogram (uniform weight if empty list)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here

bsel = b != 0

# compute dual constraints violation
Viol = alpha0[:, None] + beta0[None, :] - M
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can we change the name of the variable here by something more explicit, e.g. constraints_violation ?

# compute dual constraints violation
Viol = alpha0[:, None] + beta0[None, :] - M

# Compute worst violation per line and columns
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

worst -> largest

test/test_ot.py Outdated
np.testing.assert_almost_equal(cost1, log['cost'])
check_duality_gap(a, b, M, G, log['u'], log['v'], log['cost'])

viol = log['u'][:, None] + log['v'][None, :] - M
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here


def emd(a, b, M, numItermax=100000, log=False, dense=True):
def center_ot_dual(alpha0, beta0, a=None, b=None):
r"""Center dual OT potentials wrt theirs weights
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

wrt -> wrt.

r"""Center dual OT potentials wrt theirs weights
The main idea of this function is to find unique dual potentials
that ensure some kind of centering/fairness. It will help have
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you be more explicit on the meaning of 'centering/fairness' ?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

have -> having

@rflamary
Copy link
Collaborator Author

Thank you @ncourty i took care of your comments. We now have a warning in the emd_c function and the reason why we need to find proper dual potential is discussed.

@rflamary rflamary merged commit 81e9d42 into master Feb 3, 2020
@rflamary rflamary changed the title [WIP] Bug on dual potential (constraint violation when some weights are 0) [MRG] Bug on dual potential (constraint violation when some weights are 0) Apr 20, 2020
@rflamary rflamary deleted the bug_dual branch April 21, 2020 16:27
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants